% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

load('no_policy')

%% Generate random shocks

options = optimset;
T = 5001;               % simulation periods

%rng(1)                   % same random shock  
%shocks = normrnd(0,1,T+1,1);  % random shock
%Exi = zeros(T+1);  % stochastic steady state

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);

% obtain simulated rbv from 2-function solution
%simulated_rbv = dlmread(fullfile('from_matlab_to_fortran/simulated_rbv.txt'));

%% Simulation, stochastic steady state (shock=0)

shocks = zeros(T+1,1);

kv(1)  = k_ss;
rbv(1) = rb_ss;
nv(1)  = n_ss;
lv(1)  = kv(1)/nv(1);
av(1)  = a_ss;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);
        
        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));
        
        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            
            if pibv(t) > 0
                rhseev(t) = exp(x*solved_eta1_rhsee);
                rbv(t)    = exp(x*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x*solved_eta3_rhsee);
                rbv(t)    = exp(x*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end
            
        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x*solved_eta2_rhsee);
            rbv(t)    = exp(x*solved_eta2_rbar);
            runv(t) = 1;
        end
        
        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);
        
        erk_next = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                 + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

%        rbv(t) = fzero('func_Rb',simulated_rbv(t),options,lv(t),erk_next,lambda,gamma,delta,sigma_rk);   % deposit interest rate
        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk_next,lambda,gamma,delta,sigma_rk);   % deposit interest rate
        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        prv(t)      = normcdf((rkhats_next - erk_next) / sigma_rk);

%        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
%        z_temp     = (rkhats_next - erk_next) / sigma_rk;
%        cdf = normcdf(z_temp);
%        prv(t) = cdf;

    end
    % End of T period simulation
    
stss_n  = nv(5000);
stss_k  = kv(5000);
stss_rb = rbv(5000);

stss_y  = yv(5000);
stss_c  = cv(5000);
stss_i  = iv(5000);
stss_h  = hv(5000);
stss_l  = lv(5000);
stss_pib  = pibv(5000);
stss_pr   = prv(5000);
stss_rkhat   = rkhat(5000);
stss_rkhats  = rkhats(5000);

return
